/************************************************************************
File uts_ls_intf_openmp - utility routines for interactions with linear solver
                   implemented for openmp

Contains definitions of routines:


  utr_create_assemble_stiff_mat - to create element stiffness matrices
                                 and assemble them to the global SM
              (it calls implementations for particular platfroms)

  utr_create_assemble_stiff_mat_openmp - to create element stiffness matrices
                                 and assemble them to the global SM

------------------------------
History:
	08.2008 - Krzysztof Banas, pobanas@cyf-kr.edu.pl, initial version
*************************************************************************/

#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<math.h>
#include<assert.h>
#include<signal.h>

#include<omp.h>

/* interface for all mesh manipulation modules */
#include "mmh_intf.h"

/* interface for all approximation modules */
#include "aph_intf.h"

#ifdef PARALLEL
/* interface for parallel mesh manipulation modules */
#include "mmph_intf.h"	

/* interface with parallel communication library */
#include "pch_intf.h"
#endif

/* interface for general purpose utilities - for all problem dependent modules*/
#include "uth_intf.h"

/* interface for linear algebra packages */
#include "lin_alg_intf.h"

/* interface for control parameters and some general purpose functions */
/* from problem dependent module */
#include "pdh_control_intf.h"

#include "pdh_intf.h"

// maximal number of solution components
#define UTC_MAXEQ PDC_MAXEQ


/*------------------------------------------------------------
 utr_create_assemble_stiff_mat_openmp - to create element stiffness matrices
                                 and assemble them to the global SM
------------------------------------------------------------*/
int utr_create_assemble_stiff_mat_openmp(
  int Problem_id, 
  int Level_id, 
  int Comp_type,         /* in: indicator for the scope of computations: */
  //extern const int PDC_NO_COMP  ; /* do not compute stiff matrix and rhs vector */
  //extern const int PDC_COMP_SM  ; /* compute entries to stiff matrix only */
  //extern const int PDC_COMP_RHS ; /* compute entries to rhs vector only */
  //extern const int PDC_COMP_BOTH; /* compute entries for sm and rhsv */
  int Mth_control, // parameter to control multithreading 
  // PDC_MTH_SEQUENTIAL = 0
  // PDC_MTH_OPENMP = 1
  // PDC_MTH_OPENCL_GPU = 2
  // PDC_MTH_OPENCL_CPU = 3 ?
  // PDC_MTH_OPENCL_PHI = 4 ?
  int Nr_int_ent,
  int* L_int_ent_type,
  int* L_int_ent_id,
  int Max_dofs_int_ent,
  int* L_dof_elem_to_struct, // list of solver IDs for elements
  int* L_dof_face_to_struct, // list of solver IDs for faces
  int* L_dof_edge_to_struct, // list of solver IDs for edges
  int* L_dof_vert_to_struct  // list of solver IDs for vertices
					 );


/*------------------------------------------------------------
 utr_create_assemble_stiff_mat - to create element stiffness matrices
                                 and assemble them to the global SM
------------------------------------------------------------*/
int utr_create_assemble_stiff_mat(
  int Problem_id, 
  int Level_id, 
  int Comp_type,         /* in: indicator for the scope of computations: */
  //extern const int PDC_NO_COMP  ; /* do not compute stiff matrix and rhs vector */
  //extern const int PDC_COMP_SM  ; /* compute entries to stiff matrix only */
  //extern const int PDC_COMP_RHS ; /* compute entries to rhs vector only */
  //extern const int PDC_COMP_BOTH; /* compute entries for sm and rhsv */
  int Nr_int_ent,
  int* L_int_ent_type,
  int* L_int_ent_id,
  int Max_dofs_int_ent,
  int* L_dof_elem_to_struct, // list of solver IDs for elements
  int* L_dof_face_to_struct, // list of solver IDs for faces
  int* L_dof_edge_to_struct, // list of solver IDs for edges
  int* L_dof_vert_to_struct  // list of solver IDs for vertices
)
{

  int mth_control = PDC_MTH_OPENMP;
  utr_create_assemble_stiff_mat_openmp(Problem_id, Level_id, Comp_type, mth_control,
					 Nr_int_ent, L_int_ent_type,
					 L_int_ent_id, Max_dofs_int_ent,
					 L_dof_elem_to_struct,
					 L_dof_face_to_struct,
					 L_dof_edge_to_struct,
					 L_dof_vert_to_struct);
    

  return(1);
}


/*------------------------------------------------------------
 utr_create_assemble_stiff_mat_openmp - to create element stiffness matrices
                                 and assemble them to the global SM
------------------------------------------------------------*/
int utr_create_assemble_stiff_mat_openmp(
  int Problem_id, 
  int Level_id, 
  int Comp_type,         /* in: indicator for the scope of computations: */
  //extern const int PDC_NO_COMP  ; /* do not compute stiff matrix and rhs vector */
  //extern const int PDC_COMP_SM  ; /* compute entries to stiff matrix only */
  //extern const int PDC_COMP_RHS ; /* compute entries to rhs vector only */
  //extern const int PDC_COMP_BOTH; /* compute entries for sm and rhsv */
  int Mth_control, // parameter to control multithreading 
  // PDC_MTH_SEQUENTIAL = 0
  // PDC_MTH_OPENMP = 1
  // PDC_MTH_OPENCL_GPU = 2
  // PDC_MTH_OPENCL_CPU = 3 ?
  // PDC_MTH_OPENCL_PHI = 4 ?
  int Nr_int_ent,
  int* L_int_ent_type,
  int* L_int_ent_id,
  int Max_dofs_int_ent,
  int* L_dof_elem_to_struct, // list of solver IDs for elements
  int* L_dof_face_to_struct, // list of solver IDs for faces
  int* L_dof_edge_to_struct, // list of solver IDs for edges
  int* L_dof_vert_to_struct  // list of solver IDs for vertices
)
{


/*++++++++++++++++ executable statements ++++++++++++++++*/

/*kbw*/
  printf("In utr_create_assemble_stiff_mat_openmp: Problem_id %d, Level_id %d\n",
	 Problem_id, Level_id);
  printf("Comp_type %d, Mth_control %d, Nr_int_ent %d, Max_dofs_int_ent %d\n",
	 Comp_type, Mth_control, Nr_int_ent, Max_dofs_int_ent);
/*kew*/
  
  // divide list into elements and faces
  int nr_elems=0;  
  int int_entity;
  int previous_type = PDC_ELEMENT;
  int ok = 1;
  for(int_entity=0;int_entity<Nr_int_ent;int_entity++){
    if(L_int_ent_type[int_entity]==PDC_ELEMENT){
      nr_elems++;
      if(previous_type != PDC_ELEMENT) ok = 0;
    }
    else previous_type = PDC_FACE;
  }
  int nr_faces =  Nr_int_ent - nr_elems;

/*kbw*/
  printf("nr_elems %d, nr_faces %d, ok %d\n", nr_elems, nr_faces, ok);
/*kew*/

  if(ok!=1){
    printf("Elements are not first on the list to integrate. \nOpenMP is not optimized for that case.\n");
    //exit(-1);
  }

#pragma omp parallel if(Mth_control==PDC_MTH_OPENMP) default(none)	\
  firstprivate(Problem_id, Level_id, Comp_type,  nr_elems, Nr_int_ent, L_int_ent_type, \
	       L_int_ent_id, Max_dofs_int_ent, L_dof_elem_to_struct,	\
	       L_dof_face_to_struct, L_dof_edge_to_struct,  L_dof_vert_to_struct )
  {
    
    int nrdofs_glob, max_nrdofs, nr_dof_ent, nr_levels, posglob, nrdofs_int_ent;
    int l_dof_ent_id[PDC_MAX_DOF_PER_INT], l_dof_ent_nrdof[PDC_MAX_DOF_PER_INT];
    int l_dof_ent_posglob[PDC_MAX_DOF_PER_INT];
    int l_dof_ent_type[PDC_MAX_DOF_PER_INT];
    double *x_ini, normb;
    int i,j,k, iaux, kaux, intent, ibl, ient, nrdofbl, ini_zero;
    int level_id=0;
    char rewrite;
      
    /* allocate memory for the stiffness matrices and RHS */
    max_nrdofs = Max_dofs_int_ent;
    double *stiff_mat = (double *)malloc(max_nrdofs*max_nrdofs*sizeof(double));
    double *rhs_vect = (double *)malloc(max_nrdofs*sizeof(double));
    //double stiff_mat[100];
    //double rhs_vect[10];

    memset(stiff_mat, 0, max_nrdofs*max_nrdofs*sizeof(double));
    memset(rhs_vect, 0, max_nrdofs*sizeof(double));

/*kbw
#pragma omp critical(printing)
{
  printf("In utr_create_assemble_stiff_mat before loop over entities\n");
  printf("thread_id %d, num_threads %d\n",
	 omp_get_thread_num(),	 omp_get_num_threads() );
}
/*kew*/

    int flag = 0;

#pragma omp for // schedule(static)
    /* compute local stiffness matrices - for elements !!! */
    for(intent=0;intent<nr_elems;intent++){

      //for(intent=0;intent<0;intent++){
      //for(intent=nr_elems;intent<Nr_int_ent;intent++){
      //for(intent=0;intent<Nr_int_ent;intent++){
      
      int nrdfobl;
      int idofent;
      int nr_dof_ent = PDC_MAX_DOF_PER_INT;
      int nrdofs_int_ent = max_nrdofs;
      int l_bl_id[PDC_MAX_DOF_PER_INT], l_bl_nrdofs[PDC_MAX_DOF_PER_INT];


/*kbw*/
      if(flag==0){
#pragma omp critical(printing)
 {
  printf("In utr_create_assemble_stiff_mat before calling pdr_comp_stiff_mat\n");
  printf("intent %d, last %d, type %d, id %d, thread_id %d, num_threads %d\n",
	 intent, Nr_int_ent-1, L_int_ent_type[intent], L_int_ent_id[intent], 
	 omp_get_thread_num(), omp_get_num_threads() );
 }
      }
      flag = 1;
/*kew*/
      
  
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// ****************************************************************
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//#pragma omp critical(pdr_comp_stiff_mat)
      {
      pdr_comp_stiff_mat(Problem_id, L_int_ent_type[intent], 
			 L_int_ent_id[intent], Comp_type, NULL,
			 &nr_dof_ent,l_dof_ent_type,l_dof_ent_id,l_dof_ent_nrdof, 
			 &nrdofs_int_ent, stiff_mat, rhs_vect, &rewrite);
    /*kbw
#pragma omp critical(printing)
      {
	if(l_bl_id[1]==1){
	  //if(L_int_ent_type[intent]==PDC_ELEMENT) {
      printf("utr_create_assemble: before assemble: Solver_id %d, level_id %d, sol_typ %d\n", 
	     Problem_id, level_id, Comp_type);
      int ibl,jbl,pli,plj,nri,nrj,nrdof,jaux;
      printf("ient %d, int_ent_id %d, nr_dof_ent %d\n", 
	     intent, L_int_ent_id[intent], nrdofbl);
      pli = 0; nrdof=0;
      for(ibl=0;ibl<nrdofbl; ibl++) nrdof+=l_bl_nrdofs[ibl];
      for(ibl=0;ibl<nrdofbl; ibl++){
	printf("bl_id %d, bl_nrdof %d\n",
	  l_bl_id[ibl],l_bl_nrdofs[ibl]);
	nri = l_bl_nrdofs[ibl];
	plj=0;
	for(jbl=0;jbl<nrdofbl;jbl++){
	  printf("Stiff_mat (blocks %d:%d)\n",jbl,ibl);
	  nrj = l_bl_nrdofs[jbl];
	  for(i=0;i<nri;i++){
   	    jaux = plj+(pli+i)*nrdof;
	    for(j=0;j<nrj;j++){
	      printf("%20.15lf",stiff_mat[jaux+j]);
	      if(stiff_mat[jaux+j]< -1.e20){
		getchar(); getchar(); getchar();
	      }
	    }
	    printf("\n");
	  }
	  plj += nrj;
	}
	printf("Rhs_vect:\n");
	for(i=0;i<nri;i++){
	  printf("%20.15lf",rhs_vect[pli+i]);
	}
	printf("\n");
	pli += nri;    
      }
      getchar();
    }
      }
/*kew*/
      }

      nrdofbl = nr_dof_ent;
      
#ifdef DEBUG_SIM
      if(nrdofs_int_ent>max_nrdofs){
	printf("Too small arrays stiff_mat and rhs_vect passed to comp_el_stiff_mat\n");
	printf("from sir_create in sis_mkb. %d < %d. Exiting !!!",
	       max_nrdofs, nrdofs_int_ent);
	exit(-1);
      }
#endif
      
      
      for(idofent=0;idofent<nr_dof_ent;idofent++){
	int dof_ent_id = l_dof_ent_id[idofent];
	int dof_ent_type = l_dof_ent_type[idofent];
	int dof_ent_nrdofs = l_dof_ent_nrdof[idofent];
	int dof_struct_id;
	
	if(dof_ent_type == PDC_ELEMENT ) {
	  
	  dof_struct_id = L_dof_elem_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_elem_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	} else if(dof_ent_type == PDC_FACE ) {
	  
	  dof_struct_id = L_dof_face_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_face_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	} else if(dof_ent_type == PDC_EDGE ) {
	  
	  dof_struct_id = L_dof_edge_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_edge_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	} else if(dof_ent_type == PDC_VERTEX ) {
	  
	  dof_struct_id = L_dof_vert_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_vert_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	}
	
	/* blocks within solver are offset 1 */
	l_bl_id[idofent] = dof_struct_id + 1;
	l_bl_nrdofs[idofent] = dof_ent_nrdofs;
      }
      
    /*kbw
#pragma omp critical(printing)
      {
	if(l_bl_id[1]==1){
	  //if(L_int_ent_type[intent]==PDC_ELEMENT) {
      printf("utr_create_assemble: before assemble: Solver_id %d, level_id %d, sol_typ %d\n", 
	     Problem_id, level_id, Comp_type);
      int ibl,jbl,pli,plj,nri,nrj,nrdof,jaux;
      printf("ient %d, int_ent_id %d, nr_dof_ent %d\n", 
	     intent, L_int_ent_id[intent], nrdofbl);
      pli = 0; nrdof=0;
      for(ibl=0;ibl<nrdofbl; ibl++) nrdof+=l_bl_nrdofs[ibl];
      for(ibl=0;ibl<nrdofbl; ibl++){
	printf("bl_id %d, bl_nrdof %d\n",
	  l_bl_id[ibl],l_bl_nrdofs[ibl]);
	nri = l_bl_nrdofs[ibl];
	plj=0;
	for(jbl=0;jbl<nrdofbl;jbl++){
	  printf("Stiff_mat (blocks %d:%d)\n",jbl,ibl);
	  nrj = l_bl_nrdofs[jbl];
	  for(i=0;i<nri;i++){
   	    jaux = plj+(pli+i)*nrdof;
	    for(j=0;j<nrj;j++){
	      printf("%20.15lf",stiff_mat[jaux+j]);
	      if(stiff_mat[jaux+j]< -1.e20){
		getchar(); getchar(); getchar();
	      }
	    }
	    printf("\n");
	  }
	  plj += nrj;
	}
	printf("Rhs_vect:\n");
	for(i=0;i<nri;i++){
	  printf("%20.15lf",rhs_vect[pli+i]);
	}
	printf("\n");
	pli += nri;    
      }
      //getchar();
      }
      }
/*kew*/
    
#pragma omp critical(assembling)
	{
	  pdr_assemble_local_stiff_mat(Problem_id, level_id, Comp_type,
				       nrdofbl, l_bl_id, l_bl_nrdofs, 
				       stiff_mat, rhs_vect, &rewrite);
	}
	
    	
      } /* end loop over integration entities: ient */
    
/* #pragma omp barrier */
    flag = 0;

#pragma omp for // schedule(static)
    /* compute local stiffness matrices - for faces !!! */
    for(intent=nr_elems;intent<Nr_int_ent;intent++){
      
      //for(intent=0;intent<0;intent++){
      //for(intent=0;intent<nr_elems;intent++){
      //for(intent=0;intent<Nr_int_ent;intent++){
      
      int nrdfobl;
      int idofent;
      int nr_dof_ent = PDC_MAX_DOF_PER_INT;
      int nrdofs_int_ent = max_nrdofs;
      int l_bl_id[PDC_MAX_DOF_PER_INT], l_bl_nrdofs[PDC_MAX_DOF_PER_INT];


/*kbw*/
      if(flag==0){
  printf("In utr_create_assemble_stiff_mat before calling pdr_comp_stiff_mat\n");
  printf("intent %d, type %d, id %d, thread_id %d, num_threads %d\n",
	 intent, L_int_ent_type[intent], L_int_ent_id[intent], omp_get_thread_num(),
	 omp_get_num_threads() );
      }
      flag = 1;
/*kew*/
      
  
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// ****************************************************************
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//#pragma omp critical(pdr_comp_stiff_mat)
      {
      pdr_comp_stiff_mat(Problem_id, L_int_ent_type[intent], 
			 L_int_ent_id[intent], Comp_type, NULL,
			 &nr_dof_ent,l_dof_ent_type,l_dof_ent_id,l_dof_ent_nrdof, 
			 &nrdofs_int_ent, stiff_mat, rhs_vect, &rewrite);
      }

      nrdofbl = nr_dof_ent;
      
#ifdef DEBUG_SIM
      if(nrdofs_int_ent>max_nrdofs){
	printf("Too small arrays stiff_mat and rhs_vect passed to comp_el_stiff_mat\n");
	printf("from sir_create in sis_mkb. %d < %d. Exiting !!!",
	       max_nrdofs, nrdofs_int_ent);
	exit(-1);
      }
#endif
      
      
      for(idofent=0;idofent<nr_dof_ent;idofent++){
	int dof_ent_id = l_dof_ent_id[idofent];
	int dof_ent_type = l_dof_ent_type[idofent];
	int dof_ent_nrdofs = l_dof_ent_nrdof[idofent];
	int dof_struct_id;
	
	if(dof_ent_type == PDC_ELEMENT ) {
	  
	  dof_struct_id = L_dof_elem_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_elem_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	} else if(dof_ent_type == PDC_FACE ) {
	  
	  dof_struct_id = L_dof_face_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_face_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	} else if(dof_ent_type == PDC_EDGE ) {
	  
	  dof_struct_id = L_dof_edge_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_edge_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	} else if(dof_ent_type == PDC_VERTEX ) {
	  
	  dof_struct_id = L_dof_vert_to_struct[dof_ent_id];
	  
#ifdef DEBUG
	  if( L_dof_vert_to_struct[dof_ent_id] == -1){
	    printf("Error 347294 in sir_create!!! Exiting\n");
	    exit(-1);
	  }
#endif
	  
	}
	
	/* blocks within solver are offset 1 */
	l_bl_id[idofent] = dof_struct_id + 1;
	l_bl_nrdofs[idofent] = dof_ent_nrdofs;
      }
      
    /*kbw
#pragma omp critical(printing)
      {
	if(l_bl_id[1]==1){
    //if(L_int_ent_id[intent]!=-1) {
    int i=6;
    int solver_id = pdr_ctrl_i_params(Problem_id,i);
      printf("utr_create_assemble: before assemble: Solver_id %d, level_id %d, sol_typ %d\n", 
	     Problem_id, level_id, Comp_type);
      int ibl,jbl,pli,plj,nri,nrj,nrdof,jaux;
      printf("ient %d, int_ent_id %d, nr_dof_ent %d\n", 
	     intent, L_int_ent_id[intent], nrdofbl);
      pli = 0; nrdof=0;
      for(ibl=0;ibl<nrdofbl; ibl++) nrdof+=l_bl_nrdofs[ibl];
      for(ibl=0;ibl<nrdofbl; ibl++){
	printf("bl_id %d, bl_nrdof %d\n",
	  l_bl_id[ibl],l_bl_nrdofs[ibl]);
	nri = l_bl_nrdofs[ibl];
	plj=0;
	for(jbl=0;jbl<nrdofbl;jbl++){
	  printf("Stiff_mat (blocks %d:%d)\n",jbl,ibl);
	  nrj = l_bl_nrdofs[jbl];
	  for(i=0;i<nri;i++){
   	    jaux = plj+(pli+i)*nrdof;
	    for(j=0;j<nrj;j++){
	      printf("%20.15lf",stiff_mat[jaux+j]);

	      if(stiff_mat[jaux+j]< -1.e20){
		getchar(); getchar(); getchar();
	      }

	    }
	    printf("\n");
	  }
	  plj += nrj;
	}
	printf("Rhs_vect:\n");
	for(i=0;i<nri;i++){
	  printf("%20.15lf",rhs_vect[pli+i]);
	}
	printf("\n");
	pli += nri;    
      }
      getchar();
    }
      }
/*kew*/
    
#pragma omp critical(assembling)
	{
	  pdr_assemble_local_stiff_mat(Problem_id, level_id, Comp_type,
				       nrdofbl, l_bl_id, l_bl_nrdofs, 
				       stiff_mat, rhs_vect, &rewrite);
	}
	
    	
      } /* end loop over integration entities: ient */
    
    free(stiff_mat);
    free(rhs_vect);
  
  } // the end of parallel region

  return(1);
}


